 clear
    close all
     fid = fopen('cut1dummy.inp','r');

    keyword = 'none';
    elem_count = 1; %counts number of elements in structure
    loadcount = 1;
    nlgeom = 'NA';

    
    while ~feof(fid)
        line = upper(fgetl(fid)); % Make all lines upper case for ease
        
        if startsWith(line,'**') || isempty(line) % Process Comment Line
            continue; % Ignore comments
            
        elseif startsWith(line,'*') % Process Keyword Line
                
            if startsWith(line, '*NODE') == 1 && contains(line, 'PRINT') == 0 && contains(line, 'OUTPUT') == 0
                    keyword = 'NODE';
                elseif startsWith(line, '*ELEMENT') == 1
                    keyword = 'ELEMENT';
                    el_info = strsplit(line, ',');
                    type = extractAfter(el_info(2), '=');
                elseif startsWith(line, '*SOLID SECTION') == 1
                    keyword = 'SOLID SECTION';
                    ss_info = strsplit(line, ',');
                elseif startsWith(line, '*MATERIAL') == 1
                    keyword = '*MATERIAL';
                elseif startsWith(line, '*ELASTIC') == 1
                    keyword = 'ELASTIC';
                elseif startsWith(line, '*STEP') == 1 && contains(line, 'END') == 0
                    keyword = 'STEP';
                     
                     step_info = strsplit(line, ',');
                     if contains(step_info, 'NLGEOM')
                         nlgeom = extractAfter(step_info(3), '=');
                     end
                        if length(step_info) > 1
                            name = extractAfter(step_info(2), '=');
%                             nlgeom = extractAfter(step_info(3), '=');
                        end                   

                elseif startsWith(line, '*STATIC') == 1
                    keyword = 'STATIC';
                    static_trig = 1;
                elseif startsWith(line, '*BOUNDARY') == 1
                    keyword = 'BOUNDARY';
                elseif startsWith(line, '*CLOAD') == 1
                    keyword = 'CLOAD';
                elseif startsWith(line, '*DLOAD') == 1
                    keyword = 'DLOAD';
                elseif  startsWith(line, '*ELSET') == 1
                    keyword = 'ELSET';
            else
                keyword = 'none';
            end

        else % Process Data Line (will differ depending on keyword)
            switch keyword
                case 'NODE'
                    % Record node coordinates
                    Nint = str2num(line);
                    Model.coordinates(Nint(1), 1) = Nint(2); %stores coordinates into large matrix [#nodes,2]
                    Model.coordinates(Nint(1), 2) = Nint(3);
                    Model.coordinates(Nint(1), 3) = Nint(4);
                    
                case 'ELEMENT'
                    % Record new element definitions
                    type = extractAfter(el_info(2), '=');
                    
                    E_int = str2num(line);
                    elem_num = length(E_int); %counts number of nodes per element
                    
                    for i = 1:elem_num-1
                       Model.connectivity(elem_count, i) = E_int(i+1);
                    end
                    
                    elem_count = elem_count + 1;

            end
        end
    end

    
    fileID = fopen('uvs_1cut.txt','w');
    %fileID2 = fopen('flux_density.txt','w');
 
    xtot = max(Model.coordinates(:,1)) - min(Model.coordinates(:,1));
    ytot = max(Model.coordinates(:,2)) - min(Model.coordinates(:,2));
    ztot = max(Model.coordinates(:,3)) - min(Model.coordinates(:,3));
    ctot = [xtot/2, ytot/2, ztot/2];
    for i = 1:size(Model.connectivity,1)

        n1 = Model.connectivity(i,1); %element nodes that I need
        n2 = Model.connectivity(i,2);
        n3 = Model.connectivity(i,3);
        n4 = Model.connectivity(i,4);
        n5 = Model.connectivity(i,5);
        n6 = Model.connectivity(i,6);
        n7 = Model.connectivity(i,7);
        n8 = Model.connectivity(i,8);
        
        x1 = Model.coordinates(n1,1); %Grab coordinates per node
        y1 = Model.coordinates(n1,2);
        z1 = Model.coordinates(n1,3);
        x2 = Model.coordinates(n2,1);
        y2 = Model.coordinates(n2,2);
        z2 = Model.coordinates(n2,3);
        x3 = Model.coordinates(n3,1);
        y3 = Model.coordinates(n3,2);
        z3 = Model.coordinates(n3,3);
        x4 = Model.coordinates(n4,1);
        y4 = Model.coordinates(n4,2);
        z4 = Model.coordinates(n4,3);
        x5 = Model.coordinates(n5,1);
        y5 = Model.coordinates(n5,2);
        z5 = Model.coordinates(n5,3);
        x6 = Model.coordinates(n6,1);
        y6 = Model.coordinates(n6,2);
        z6 = Model.coordinates(n6,3);
        x7 = Model.coordinates(n7,1);
        y7 = Model.coordinates(n7,2);
        z7 = Model.coordinates(n7,3);
        x8 = Model.coordinates(n8,1);
        y8 = Model.coordinates(n8,2);
        z8 = Model.coordinates(n8,3);
        
        avx = mean([x1,x2,x3,x4,x5,x6,x7,x8]);
        avy = mean([y1,y2,y3,y4,y5,y6,y7,y8]);
        avz = mean([z1,z2,z3,z4,z5,z6,z7,z8]);
        
        cent = [avx,avy,avz];
        
        Model.centroids(i,:) = cent;
        
    end
   
%fid = fopen('1.5875radmagnet_extrafinemesh.txt','r'); %magnet radius = 1.5875, h = 3.175, extra fine mesh
fid = fopen('arrows0628.txt','r'); %(normal_standard)
%fid = fopen('PMArrows_small.txt','r'); (small mag)

%fid = fopen('3magArrows.txt','r');

i = 1;
while ~feof(fid)
    line_i = upper(fgetl(fid)); % Make all lines upper case for ease
    if startsWith(line_i,'%') || isempty(line) % Process Comment Line
            continue; % Ignore comments
    else
        line = str2num(line_i);
        %misalignment = 0.0;
        if line(1) > -11 && line(1) < 11 && line(2) > -11 && line(2) < 11
        x(i) = line(1);
        y(i) = line(2);
        z(i) = line(3);
        Mx(i) = line(4);
        My(i) = line(5);
        Mz(i) = line(6);
        vec = [Mx(i), My(i), Mz(i)];
        uvec(i,:) = vec / sqrt(vec(1)^2 + vec(2)^2 + vec(3)^2);
        Bx(i) = uvec(i,1);
        By(i) = uvec(i,2);
        Bz(i) = uvec(i,3);
        i = i+1;
        else
            continue
        end
    end
end
B(:,1) = Bx;
B(:,2) = By;
B(:,3) = Bz;
loc(:,1) = x;
loc(:,2) = y;
loc(:,3) = z;


for m = 1:size(Model.centroids,1)
    xx(m) = Model.centroids(m,1);
    yy(m) = Model.centroids(m,2);
    zz(m) = Model.centroids(m,3);
    Idx = knnsearch(loc,Model.centroids(m,:));
    vec = B(Idx,:);
%     
%     if vec(3) < 0.000000
%         vec = [0 0 0];
%     else
        vec = vec/sqrt(vec(1)^2 + vec(2)^2 + vec(3)^2); 
%     end

    unitvec = vec;
    
   
    fprintf(fileID,'%f %f %f \n', vec);

    pBx(m) = unitvec(1);
    pBy(m) = unitvec(2);
    pBz(m) = unitvec(3);
end

quiver3(xx,yy,zz,pBx,pBy,pBz,'linewidth',2)
% hold on
% theta = 0:pi/50:2*pi;
% x_c = 6.35*cos(theta);
% y_c = 6.35*sin(theta);
% %plot(x_c,y_c)
% fill(x_c,y_c,[0.7 0.7 0.7],'FaceAlpha',0.3,'EdgeAlpha',0.0)
daspect([1 1 1])
grid off



